# Replication of analysis for
# Monitoring Coalition Partners in the EU: Strategic Committee Appointments in the European Parliament
# Pit Rieger --- November 2025
library(here) # v1.0.1
library(tidyverse) # v2.0.0
library(modelsummary) #2.2.0
library(marginaleffects) # v0.22.0
library(fixest) # v0.13.2
library(latex2exp) # v0.9.6
library(ggdist) # v3.3.2
library(distributional) # v0.5.0
library(stargazer) # v5.2.3
library(interflex) # v1.2.6
library(dotwhisker) # v0.8.2
library(survival) # v3.6-4
library(tinytable) # v0.4.0
sessionInfo()

# read data
load(here("Rieger2025_data.RData"))

# Main analysis ====
  ## Committee appointment ----
    ### Fit models ----
      # Model 1
      mCOM.policyFE1 = feglm(
        member ~
          oppoutin + oppoutin:delta_policy
        |
          1,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyFE1)

      #### Model 2
      mCOM.policyFEcabparty = feglm(
        member ~
          oppoutin + oppoutin:delta_policy
        |
          cabparty,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyFEcabparty)

      # Model 3
      mCOM.policyFEcabpartycom = feglm(
        member ~
          oppoutin + oppoutin:delta_policy
        |
          cabparty + title_short,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyFEcabpartycom)

      # Model 4
      mCOM.policyFEcabpartycom_term = feglm(
        member ~
          oppoutin + oppoutin:delta_policy
        |
          cabparty^EPterm + title_short^EPterm,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyFEcabpartycom_term)

    ### Generate Table 1 ----
    models_policy_COM = list(mCOM.policyFE1,
                             mCOM.policyFEcabparty,
                             mCOM.policyFEcabpartycom,
                             mCOM.policyFEcabpartycom_term)
    FE_add = data.frame(matrix(c("Party $\\times$ Cabinet FE","", "\\checkmark", "\\checkmark", "\\checkmark",
                                 "Committee FE", "", "", "\\checkmark", "\\checkmark",
                                 "Committee $\\times$ term FE", "", "", "", "\\checkmark",
                                 "Party $\\times$ cabinet $\\times$ term FE", "", "", "", "\\checkmark",
                                 "Log. Lik.", round(unlist(lapply(models_policy_COM, logLik)), 3)),
                               ncol = length(models_policy_COM) + 1, byrow = T))
    attr(FE_add, "position") = 12 + 1:nrow(FE_add)

    regtab = modelsummary(models_policy_COM,
                          #output = "latex",
                          coef_map = cmap,
                          gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                          add_rows = FE_add,
                          escape=FALSE,
                          stars = T,
                          title = "Binary logistic regression estimates. DV: Committee appointment. Standard errors in parentheses.\\label{tab:regtab_deltapolicy_COM}")
    regtab
    save_tt(regtab, here("tab1.txt"), overwrite = T)
    rm(regtab, FE_add, models_policy_COM)

    ### Generate Figure 1 ----
    pdat = plot_slopes(mCOM.policyFEcabpartycom_term, variables = "delta_policy", by = "oppoutin", draw = F)
    pdat
    pdat$oppoutin = case_match(pdat$oppoutin,
                               "outsider" ~ "Outsider",
                               "opposition" ~ "Opposition",
                               "insider" ~ "Insider")
    ggplot(pdat, aes(x = oppoutin, ydist = dist_normal(estimate, std.error))) +
      geom_hline(yintercept = 0, linetype = 2) +
      stat_eye(fill = "grey80", .width = c(.90, .95), alpha = 0.5) +
      stat_eye(fill = NA, .width = c(.90, .95)) +
      theme_bw() +
      labs(x = "Party status", y = TeX("Effect of misalignment ($\\Delta_{policy}$) on committee membership"))
    ggsave(here("fig1.pdf"), width = 6.5, height = 4)

  ## Rapporteur Assignment ----
    ### Fit models ----
    mRAP.policyagreeFE1 = feglm(
      appointed ~
        member_pct +
        oppoutin*delta_policy*billdisagree_wo
      |
        1
      ,
      family = binomial,
      vcov = "iid",
      fixef.rm = "none",
      data = d_RAP,
      fixef.keep_names = T
    )
    summary(mRAP.policyagreeFE1)

    mRAP.policyagreeFEcabparty = feglm(
      appointed ~
        member_pct +
        oppoutin*delta_policy*billdisagree_wo
      |
        cabparty
      ,
      family = binomial,
      vcov = "iid",
      fixef.rm = "none",
      data = d_RAP,
      fixef.keep_names = T
    )
    summary(mRAP.policyagreeFEcabparty)

    mRAP.policyagreeFEcabpartycom = feglm(
      appointed ~
        member_pct +
        oppoutin*delta_policy*billdisagree_wo
      |
        cabparty +
        title_short
      ,
      family = binomial,
      vcov = "iid",
      fixef.rm = "none",
      data = d_RAP,
      fixef.keep_names = T
    )
    summary(mRAP.policyagreeFEcabpartycom)

    mRAP.policyagreeFEcabpartycom_term = feglm(
      appointed ~
        member_pct +
        oppoutin*delta_policy*billdisagree_wo
      |
        cabparty^EPterm +
        title_short^EPterm
      ,
      family = binomial,
      vcov = "iid",
      fixef.rm = "none",
      data = d_RAP,
      fixef.keep_names = T
    )
    summary(mRAP.policyagreeFEcabpartycom_term)

  ### Generate Table A5 ----
  models_policyagree_RAP = list(mRAP.policyagreeFE1,
                                mRAP.policyagreeFEcabparty,
                                mRAP.policyagreeFEcabpartycom,
                                mRAP.policyagreeFEcabpartycom_term)
  FE_add = data.frame(matrix(c("Party $\\times$ Cabinet FE","", "\\checkmark", "\\checkmark", "\\checkmark",
                               "Committee FE", "", "", "\\checkmark", "\\checkmark",
                               "Committee $\\times$ term FE", "", "", "", "\\checkmark",
                               "Party $\\times$ cabinet $\\times$ term FE", "", "", "", "\\checkmark",
                               "Log. Lik.", round(unlist(lapply(models_policyagree_RAP, logLik)), 3)),
                             ncol = length(models_policyagree_RAP) + 1, byrow = T))
  attr(FE_add, "position") = 26 + 1:nrow(FE_add)

  regtab = modelsummary(models_policyagree_RAP,
                        #output = "latex",
                        coef_map = cmap,
                        gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                        add_rows = FE_add,
                        escape=FALSE,
                        stars = T,
                        width = c(0.4, rep(0.6/length(models_policyagree_RAP), length(models_policyagree_RAP))),
                        title = "Binary logistic regression estimates. DV: Rapporteur assignment. Standard errors in parentheses.\\label{tab:regtab_policyagree_RAP}")# c('*' = .05, '**' = .01, "***" = .001))
  regtab
  save_tt(regtab, here("tabA5.txt"), overwrite = T)
  rm(regtab, FE_add, models_policyagree_RAP)

  ### Generate Figure 2 ----
  pdat = plot_slopes(mRAP.policyagreeFEcabpartycom_term, variables = "delta_policy", condition = list("billdisagree_wo", "oppoutin"), draw = F)
  pdat$oppoutin = case_match(pdat$oppoutin,
                             "opposition" ~ "Opposition",
                             "outsider" ~ "Outsider",
                             "insider" ~ "Insider")
  pdat %>% filter(billdisagree_wo == 1)

  ggplot(pdat, aes(x = billdisagree_wo, y = estimate, ymin = conf.low, ymax = conf.high)) +
    geom_rug(inherit.aes = F,
             data = d_RAP %>% dplyr::select(billdisagree_wo,oppoutin) %>%
               mutate(oppoutin = case_match(oppoutin,
                                            "opposition" ~ "Opposition",
                                            "outsider" ~ "Outsider",
                                            "insider" ~ "Insider")),
             mapping = aes(x = billdisagree_wo),
             alpha = 0.3) +
    facet_wrap(~oppoutin)+
    geom_ribbon(linetype = 2, color = "grey20", fill = "grey80", alpha = 0.5) +
    geom_line() +
    theme_bw() +
    theme(strip.background = element_blank(),
          strip.text = element_text(size = 11)) +
    labs(x = "Bill disagreement",
         y = TeX("Effect of misalignment ($\\Delta_{policy}$) on rapporteur assignment"),
    )
  ggsave(here("fig2.pdf"), width = 6.5, height = 4)

# Appendix ====
  ## Descriptive Statistics ----
    ### Committee appointment Table A1----
    sumvars = c("member",
                "oppoutinoutsider",
                "oppoutininsider",
                "oppoutinopposition",
                "delta_policy",
                "delta_lr",
                "delta_policy_outsider")
    cmap[sumvars]
    sumvars[!sumvars %in% colnames(d_COM)]
    dstat_COM = stargazer(as.data.frame(d_COM[,sumvars]), summary = T, covariate.labels = cmap[sumvars])
    writeLines(dstat_COM, here("tabA1.txt"))

    ### Rapporteur assignment Table A2----
    sumvars = c("appointed",
                "oppoutinoutsider",
                "oppoutininsider",
                "oppoutinopposition",
                "member_pct",
                "delta_policy",
                "delta_lr",
                "delta_policy_outsider",
                "billdisagree_wo")
    cmap[sumvars]

    sumvars[!sumvars %in% colnames(d_RAP)]
    dstat_RAP = stargazer(as.data.frame(d_RAP[,sumvars]), summary = T, covariate.labels = cmap[sumvars])
    writeLines(dstat_RAP, here("tabA2.txt"))

    ### Generate Figure A3 ----
    ggplot(d_RAP %>%
             mutate(oppoutin = case_match(oppoutin,
                                          "opposition" ~ "Opposition",
                                          "outsider" ~ "Outsider",
                                          "insider" ~ "Insider")),
           aes(y = billdisagree_wo, x = oppoutin)) +
      geom_violin(alpha = 0.2) +
      geom_point(stat = "summary",
                 fun = mean) +
      theme_bw() +
      theme(strip.background = element_blank(),
            strip.text = element_text(size = 11)) +
      labs(y = "Bill disagreement",
           x = "Party status")
    ggsave(here("figA3.pdf"), width = 6.5, height = 4)

  ## MEP survey ----
    ### Contact frequency ----
      # prep data
      d_contact = d_MEPsurvey %>% filter(!is.na(contact_natministers)) %>%
        filter(gov_isCoalition == 1) %>%
        group_by(gov, contact_natministers) %>%
        summarize(n = n()) %>%
        group_by(gov) %>%
        mutate(pct = 100*n/sum(n)) %>%
        ungroup() %>%
        mutate(gov = factor(gov),
               var = "National minister")
      d_contact =
        bind_rows(d_contact,
                  d_MEPsurvey %>% filter(!is.na(contact_natpartyleader)) %>%
                    filter(gov_isCoalition == 1) %>%
                    group_by(gov, contact_natpartyleader) %>%
                    summarize(n = n()) %>%
                    group_by(gov) %>%
                    mutate(pct = 100*n/sum(n)) %>%
                    ungroup() %>%
                    mutate(gov = factor(gov),
                           var = "National party (leadership)"))
      d_contact =
        bind_rows(d_contact,
                  d_MEPsurvey %>% filter(!is.na(contact_copartisans)) %>%
                    filter(gov_isCoalition == 1) %>%
                    group_by(gov, contact_copartisans) %>%
                    summarize(n = n()) %>%
                    group_by(gov) %>%
                    mutate(pct = 100*n/sum(n)) %>%
                    ungroup() %>%
                    mutate(gov = factor(gov),
                           var = "National party (member)"))

      # clean labels
      d_contact = d_contact %>%
        mutate(contact = coalesce(contact_natministers, contact_natpartyleader, contact_copartisans),
               contact_labels = case_when(contact == 1 ~ "week",
                                          contact == 2 ~ "month",
                                          contact == 3 ~ "3 months",
                                          contact == 4 ~ "year",
                                          contact == 5 ~ "less often",
                                          contact == 6 ~ "never")) %>%
        arrange(contact) %>%
        mutate(contact_labels = factor(contact_labels, levels = unique(contact_labels)))
      d_contact$var = factor(d_contact$var, levels = c("National party (leadership)", "National party (member)", "National minister"))

      # Generate Figure A1
      ggplot(d_contact, aes(x = contact_labels, y = pct, fill = gov)) +
        geom_bar(stat = "identity", position = "dodge", color = "grey20") +
        scale_fill_manual(values = c("grey80", "grey50")) +
        labs(x = "Frequency of contact",
             y = "% of MEPs") +
        guides(fill = guide_legend(title = "MEP from government party", direction = "horizontal")) +
        theme_bw() +
        theme(legend.position = "bottom") +
        facet_wrap(~var, ncol = 2) +
        theme(legend.position = c(0.8,0.1),
              legend.title.position = "top",
              axis.text.x = element_text(angle = 20, hjust = 1),
              strip.text = element_text(face = "bold", hjust = 0, size = 11),
              strip.background = element_blank())
      ggsave(here("figA1.pdf"), width = 6.5, height = 4)

    ### Committee choice ----
      # prep data
      d_committeechoice = d_MEPsurvey %>% filter(!is.na(committeechoice_natparty)) %>%
        filter(gov_isCoalition == 1) %>%
        group_by(gov, committeechoice_natparty) %>%
        summarize(n = n()) %>%
        group_by(gov) %>%
        mutate(pct = 100*n/sum(n),
               gov = factor(gov),
               var = "National party")
      d_committeechoice =
        bind_rows(d_committeechoice,
                  d_MEPsurvey %>% filter(!is.na(committeechoice_EPgroup)) %>%
                    filter(gov_isCoalition == 1) %>%
                    group_by(gov, committeechoice_EPgroup) %>%
                    summarize(n = n()) %>%
                    group_by(gov) %>%
                    mutate(pct = 100*n/sum(n),
                           gov = factor(gov),
                           var = "Party group")) %>%
        mutate(importance = coalesce(committeechoice_natparty, committeechoice_EPgroup)) %>%
        select(-committeechoice_natparty, -committeechoice_EPgroup)

      # Generate Figure A2
      ggplot(d_committeechoice, aes(x = importance, y = pct, fill = gov)) +
        geom_bar(stat = "identity", position = "dodge", color = "grey20") +
        facet_wrap(~var) +
        #scale_x_continuous() +
        scale_fill_manual(values = c("grey80", "grey50")) +

        scale_x_reverse(breaks = 5:1, labels = c("Not at all\nimportant", 4:2, "Extremely\nimportant")) +
        theme_bw() +
        labs(x = "Importance for committee choice",
             y = "% of MEPs") +
        guides(fill = guide_legend(title = "MEP from government party")) +
        theme_bw() +
        theme(legend.position = "bottom",
              strip.text = element_text(face = "bold", hjust = 0, size = 11),
              strip.background = element_blank())
      ggsave(here("figA2.pdf"), width = 6.5, height = 4)

    ## Remove member_pct from ----
      ### Fit models ----
      mRAP.policyagreeFE1_exclmemberpct = feglm(
        appointed ~
          oppoutin*delta_policy*billdisagree_wo
        |
          1
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFE1_exclmemberpct)

      mRAP.policyagreeFEcabparty_exclmemberpct = feglm(
        appointed ~
          oppoutin*delta_policy*billdisagree_wo
        |
          cabparty
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFEcabparty_exclmemberpct)

      mRAP.policyagreeFEcabpartycom_exclmemberpct = feglm(
        appointed ~
          oppoutin*delta_policy*billdisagree_wo
        |
          cabparty +
          title_short
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFEcabpartycom_exclmemberpct)

      mRAP.policyagreeFEcabpartycom_term_exclmemberpct = feglm(
        appointed ~
          oppoutin*delta_policy*billdisagree_wo
        |
          cabparty^EPterm +
          title_short^EPterm
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFEcabpartycom_term_exclmemberpct)

      ### Generate table A6----
      models_policyagree_RAP_exclmemberpct = list(mRAP.policyagreeFE1_exclmemberpct,
                                                  mRAP.policyagreeFEcabparty_exclmemberpct,
                                                  mRAP.policyagreeFEcabpartycom_exclmemberpct,
                                                  mRAP.policyagreeFEcabpartycom_term_exclmemberpct)
      FE_add = data.frame(matrix(c("Party $\\times$ Cabinet FE","", "\\checkmark", "\\checkmark", "\\checkmark",
                                   "Committee FE", "", "", "\\checkmark", "\\checkmark",
                                   "Committee $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Party $\\times$ cabinet $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Log. Lik.", round(unlist(lapply(models_policyagree_RAP_exclmemberpct, logLik)), 3)),
                                 ncol = length(models_policyagree_RAP_exclmemberpct) + 1, byrow = T))
      attr(FE_add, "position") = 24 + 1:nrow(FE_add)

      regtab = modelsummary(models_policyagree_RAP_exclmemberpct,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            width = c(0.4, rep(0.6/length(models_policyagree_RAP_exclmemberpct), length(models_policyagree_RAP_exclmemberpct))),
                            title = "Binary logistic regression estimates. DV: Rapporteur assignment. Standard errors in parentheses.\\label{tab:regtab_policyagree_RAP_exclmemberpct}")# c('*' = .05, '**' = .01, "***" = .001))
      regtab
      save_tt(regtab, here("tabA6.txt"), overwrite = T)
      rm(regtab, FE_add, models_policyagree_RAP_exclmemberpct)

  ## LPM ----
    ### Committee appointment ----
      # Fit models
      mCOM.policyFE1_LPM = feols(
        member ~
          oppoutin + oppoutin:delta_policy
        |
          1,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyFE1_LPM)

      mCOM.policyFEcabparty_LPM = feols(
        member ~
          oppoutin + oppoutin:delta_policy
        |
          cabparty,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyFEcabparty_LPM)

      mCOM.policyFEcabpartycom_LPM = feols(
        member ~
          oppoutin + oppoutin:delta_policy
        |
          cabparty + title_short,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyFEcabpartycom_LPM)

      mCOM.policyFEcabpartycom_term_LPM = feols(
        member ~
          oppoutin + oppoutin:delta_policy
        |
          cabparty^EPterm + title_short^EPterm,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyFEcabpartycom_term_LPM)

      # Generate table A7
      models_policy_COM_LPM = list(mCOM.policyFE1_LPM,
                                   mCOM.policyFEcabparty_LPM,
                                   mCOM.policyFEcabpartycom_LPM,
                                   mCOM.policyFEcabpartycom_term_LPM)
      FE_add = data.frame(matrix(c("Party $\\times$ Cabinet FE","", "\\checkmark", "\\checkmark", "\\checkmark",
                                   "Committee FE", "", "", "\\checkmark", "\\checkmark",
                                   "Committee $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Party $\\times$ cabinet $\\times$ term FE", "", "", "", "\\checkmark"),
                                 ncol = length(models_policy_COM_LPM) + 1, byrow = T))
      attr(FE_add, "position") = 12 + 1:nrow(FE_add)

      regtab = modelsummary(models_policy_COM_LPM,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|Within",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            title = "Linear regression estimates. DV: Committee appointment. Standard errors in parentheses.\\label{tab:regtab_deltapolicy_COM_LPM}")# c('*' = .05, '**' = .01, "***" = .001))
      regtab
      save_tt(regtab, here("tabA7.txt"), overwrite = T)
      rm(regtab, FE_add, models_policy_COM_LPM)

      # Generate Figure A4
      pdat = plot_slopes(mCOM.policyFEcabpartycom_term_LPM, variables = "delta_policy", by = "oppoutin", draw = F)
      pdat$oppoutin = case_match(pdat$oppoutin,
                                 "outsider" ~ "Outsider",
                                 "opposition" ~ "Opposition",
                                 "insider" ~ "Insider")
      ggplot(pdat, aes(x = oppoutin, ydist = dist_normal(estimate, std.error))) +
        geom_hline(yintercept = 0, linetype = 2) +
        stat_eye(fill = "grey80", .width = c(.90, .95), alpha = 0.5) +
        stat_eye(fill = NA, .width = c(.90, .95)) +

        theme_bw() +
        labs(x = "Party status", y = TeX("Effect of misalignment ($\\Delta_{policy}$) on committee membership"))
      ggsave(here("figA4.pdf"), width = 6.5, height = 4)

    ### Rapporteur assignment ----
      # Fit models
      mRAP.policyagreeFE1_LPM = feols(
        appointed ~
          member_pct +
          oppoutin*delta_policy*billdisagree_wo
        |
          1
        ,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFE1_LPM)

      mRAP.policyagreeFEcabparty_LPM = feols(
        appointed ~
          member_pct +
          oppoutin*delta_policy*billdisagree_wo
        |
          cabparty
        ,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFEcabparty_LPM)

      mRAP.policyagreeFEcabpartycom_LPM = feols(
        appointed ~
          member_pct +
          oppoutin*delta_policy*billdisagree_wo
        |
          cabparty +
          title_short
        ,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFEcabpartycom_LPM)

      mRAP.policyagreeFEcabpartycom_term_LPM = feols(
        appointed ~
          member_pct +
          oppoutin*delta_policy*billdisagree_wo
        |
          cabparty^EPterm +
          title_short^EPterm
        ,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFEcabpartycom_term_LPM)

      # Generate table A8
      models_policyagree_RAP_LPM = list(mRAP.policyagreeFE1_LPM,
                                        mRAP.policyagreeFEcabparty_LPM,
                                        mRAP.policyagreeFEcabpartycom_LPM,
                                        mRAP.policyagreeFEcabpartycom_term_LPM)

      FE_add = data.frame(matrix(c("Party $\\times$ Cabinet FE","", "\\checkmark", "\\checkmark", "\\checkmark",
                                   "Committee FE", "", "", "\\checkmark", "\\checkmark",
                                   "Committee $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Party $\\times$ cabinet $\\times$ term FE", "", "", "", "\\checkmark"),
                                 ncol = length(models_policyagree_RAP_LPM) + 1, byrow = T))
      attr(FE_add, "position") = 26 + 1:nrow(FE_add)

      regtab = modelsummary(models_policyagree_RAP_LPM,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|Within",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            width = c(0.4, rep(0.6/length(models_policyagree_RAP_LPM), length(models_policyagree_RAP_LPM))),
                            title = "Linear regression estimates. DV: Rapporteur assignment. Standard errors in parentheses.\\label{tab:regtab_policyagree_RAP_LPM}")# c('*' = .05, '**' = .01, "***" = .001))
      regtab
      save_tt(regtab, here("tabA8.txt"), overwrite = T)
      rm(regtab, FE_add, models_policyagree_RAP_LPM)

      # Generate Figure A5
      pdat = plot_slopes(mRAP.policyagreeFEcabpartycom_term_LPM, variables = "delta_policy", condition = list("billdisagree_wo", "oppoutin"), draw = F)
      pdat$oppoutin = case_match(pdat$oppoutin,
                                 "opposition" ~ "Opposition",
                                 "outsider" ~ "Outsider",
                                 "insider" ~ "Insider")

      ggplot(pdat, aes(x = billdisagree_wo, y = estimate, ymin = conf.low, ymax = conf.high)) +
        geom_rug(inherit.aes = F,
                 data = d_RAP %>% dplyr::select(billdisagree_wo,oppoutin) %>%
                   mutate(oppoutin = case_match(oppoutin,
                                                "opposition" ~ "Opposition",
                                                "outsider" ~ "Outsider",
                                                "insider" ~ "Insider")),
                 mapping = aes(x = billdisagree_wo),
                 alpha = 0.3) +
        facet_wrap(~oppoutin)+
        geom_ribbon(linetype = 2, color = "grey20", fill = "grey80", alpha = 0.5) +
        geom_line() +
        theme_bw() +
        theme(strip.background = element_blank(),
              strip.text = element_text(size = 11)) +
        labs(x = "Bill disagreement",
             y = TeX("Effect of misalignment ($\\Delta_{policy}$) on rapporteur assignment"),
        )
      ggsave(here("figA5.pdf"), width = 6.5, height = 4)

  ## Conditional logit models ----
    ### Committee appointment ----
      # fit model
      d_COM$cabpartyterm = paste0(d_COM$cabparty, "-", d_COM$EPterm)
      mCOM.policyFEcabpartyterm_clogit = clogit(
        member ~ oppoutin + oppoutin:delta_policy +
          strata(cabpartyterm),
        data = d_COM)
      summary(mCOM.policyFEcabpartyterm_clogit)

      # Generate table A9
      models_policy_COM_clogit = list(mCOM.policyFEcabpartyterm_clogit)
      FE_add = data.frame(matrix(c("Log. Lik.", round(unlist(lapply(models_policy_COM_clogit, logLik)), 3)),
                                 ncol = length(models_policy_COM_clogit) + 1, byrow = T))
      attr(FE_add, "position") = 8 + 1:nrow(FE_add)

      regtab = modelsummary(models_policy_COM_clogit,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            title = "Conditional logistic regression estimates. Replication of model 4. DV: Committee appointment. Standard errors in parentheses.\\label{tab:regtab_policy_COM_clogit}")
      regtab
      save_tt(regtab, here("tabA9.txt"), overwrite = T)
      rm(regtab, models_policyagree_RAP_clogit, FE_add)

    ### Rapporteur assignment ----
      # fit model
      d_RAP$cabpartyterm = paste0(d_RAP$cabparty, "-", d_RAP$EPterm)
      mRAP.policyagreeFEcabpartyterm_clogit <- clogit(
        appointed ~ oppoutin*delta_policy*billdisagree_wo +
          member_pct +
          strata(cabpartyterm),
        method = "efron",
        data = d_RAP
      )
      summary(mRAP.policyagreeFEcabpartyterm_clogit)

      # Generate table A10
      models_policyagree_RAP_clogit = list(mRAP.policyagreeFEcabpartyterm_clogit)
      FE_add = data.frame(matrix(c("Log. Lik.", round(unlist(lapply(models_policyagree_RAP_clogit, logLik)), 3)),
                                 ncol = length(models_policyagree_RAP_clogit) + 1, byrow = T))
      attr(FE_add, "position") = 24 + 1:nrow(FE_add)

      regtab = modelsummary(models_policyagree_RAP_clogit,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            title = "Conditional logistic regression estimates. Replication of model A4. DV: Rapporteur assignment. Standard errors in parentheses.\\label{tab:regtab_policyagree_RAP_clogit}")
      regtab
      save_tt(regtab, here("tabA10.txt"), overwrite = T)
      rm(regtab, models_policyagree_RAP_clogit, FE_add)

  ## Multiplicative interaction robustness ----
    ### Committee appointment ----
      # prep to accommodate the necessary package requirements
      d_COM_interaction = d_COM
      d_COM_interaction = d_COM_interaction %>%
        mutate(oppoutin = factor(oppoutin),
               party_id_fac = factor(paste0(EPterm , "-", cabparty)),
               title_short_fac = factor(paste0(EPterm, "-", title_short)))

      # fit
      p2 = interflex(estimator = "binning",
                     D = "outsider",
                     X = "delta_policy",
                     Y = "member",
                     data = d_COM_interaction,
                     FE = c("party_id_fac", "title_short_fac"),
                     nbins = 2,
                     na.rm = T)
      p3 = interflex(estimator = "binning",
                     D = "outsider",
                     X = "delta_policy",
                     Y = "member",
                     data = d_COM_interaction,
                     FE = c("party_id_fac", "title_short_fac"),
                     nbins = 3,
                     na.rm = T)
      p4 = interflex(estimator = "binning",
                     D = "outsider",
                     X = "delta_policy",
                     Y = "member",
                     data = d_COM_interaction,
                     FE = c("party_id_fac", "title_short_fac"),
                     nbins = 4,
                     na.rm = T)
      set.seed(111)
      pkernel = interflex(estimator = "kernel",
                          D = "outsider",
                          X = "delta_policy",
                          Y = "member",
                          data = d_COM_interaction,
                          FE = c("party_id_fac", "title_short_fac"),
                          na.rm = T)

      # Generate Figure A6
      cowplot::plot_grid(p2$figure, p3$figure, p4$figure, pkernel$figure)
      ggsave(here("figA6.pdf"), width = 13, height = 13)

    ### Rapporteur assignment ----
      # prep to accommodate the necessary package requirements
      d_RAP_interaction = d_RAP
      d_RAP_interaction = d_RAP_interaction %>%
        mutate(oppoutin = factor(oppoutin),
               outsider = factor(outsider),
               party_id_fac = factor(paste0(cabparty, "-", EPterm)),
               title_short_fac = factor(paste0(title_short, "-", EPterm)))

      # fit
      p2 = interflex(estimator = "binning",
                     D = "delta_policy",
                     X = "billdisagree_wo",
                     Y = "appointed",
                     Z = "member_pct",
                     data = d_RAP_interaction %>% filter(oppoutin == "outsider" & !is.na(billdisagree_wo)),
                     FE = c("party_id_fac", "title_short_fac"),
                     nbins = 2,
                     na.rm = F)
      p3 = interflex(estimator = "binning",
                     D = "delta_policy",
                     X = "billdisagree_wo",
                     Y = "appointed",
                     Z = "member_pct",
                     data = d_RAP_interaction %>% filter(oppoutin == "outsider" & !is.na(billdisagree_wo)),
                     FE = c("party_id_fac", "title_short_fac"),
                     nbins = 4,
                     na.rm = F)
      p4 = interflex(estimator = "binning",
                     D = "delta_policy",
                     X = "billdisagree_wo",
                     Y = "appointed",
                     Z = "member_pct",
                     data = d_RAP_interaction %>% filter(oppoutin == "outsider" & !is.na(billdisagree_wo)),
                     FE = c("party_id_fac", "title_short_fac"),
                     nbins = 5,
                     na.rm = F)
      pkernel = interflex(estimator = "kernel",
                          D = "delta_policy",
                          X = "billdisagree_wo",
                          Y = "appointed",
                          Z = "member_pct",
                          data = d_RAP_interaction %>% filter(oppoutin == "outsider" & !is.na(billdisagree_wo)),
                          FE = c("party_id_fac", "title_short_fac"),
                          na.rm = F)

      # Generate Figure A7
      cowplot::plot_grid(p2$figure, p3$figure, p4$figure, pkernel$figure)
      ggsave(here("figA7.pdf"), width = 13, height = 13)

  ## RILE positions ----
    ### Committee appointment ----
      # Fit models
      mCOM.rileFE1 = feglm(
        member ~
          oppoutin + oppoutin:delta_lr
        |
          1,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.rileFE1)

      mCOM.rileFEparty = feglm(
        member ~
          oppoutin + oppoutin:delta_lr
        |
          party_id,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.rileFEparty)

      mCOM.rileFEpartycom = feglm(
        member ~
          oppoutin + oppoutin:delta_lr
        |
          party_id + title_short,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.rileFEpartycom)

      mCOM.rileFEpartycom_term = feglm(
        member ~
          oppoutin + oppoutin:delta_lr
        |
          party_id + title_short^EPterm,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.rileFEpartycom_term)

      # Generate table A11
      models_rile_COM = list(mCOM.rileFE1,
                             mCOM.rileFEparty,
                             mCOM.rileFEpartycom,
                             mCOM.rileFEpartycom_term)
      FE_add = data.frame(matrix(c("Party FE", "", "\\checkmark", "\\checkmark", "\\checkmark",
                                   "Committee FE", "", "", "\\checkmark", "\\checkmark",
                                   "Term FE", "", "", "", "\\checkmark",
                                   "Committee $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Log. Lik.", round(unlist(lapply(models_rile_COM, logLik)), 3)),
                                 ncol = length(models_rile_COM) + 1, byrow = T))
      attr(FE_add, "position") = 12 + 1:nrow(FE_add)

      regtab = modelsummary(models_rile_COM,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            title = "Binary logistic regression estimates. DV: Committee appointment. Standard errors in parentheses.\\label{tab:regtab_deltarile_COM}")# c('*' = .05, '**' = .01, "***" = .001))
      regtab
      save_tt(regtab, here("tabA11.txt"), overwrite = T)
      rm(regtab, FE_add, models_rile_COM)

      # Generate Figure A8
      pdat = plot_slopes(mCOM.rileFEpartycom_term, variables = "delta_lr", by = "oppoutin", draw = F)
      pdat
      pdat$oppoutin = case_match(pdat$oppoutin,
                                 "outsider" ~ "Outsider",
                                 "opposition" ~ "Opposition",
                                 "insider" ~ "Insider")
      ggplot(pdat, aes(x = oppoutin, ydist = dist_normal(estimate, std.error))) +
        geom_hline(yintercept = 0, linetype = 2) +
        stat_eye(fill = "grey60", .width = c(.90, .95), alpha = 0.5) +
        stat_eye(fill = NA, .width = c(.90, .95)) +
        theme_bw() +
        labs(x = "Party status", y = TeX("Effect of misalignment ($\\Delta_{rile}$) on committee membership"))
      ggsave(here("figA8.pdf"), width = 6.5, height = 4)

    ### Rapporteur assignment ----
      # Fit models
      mRAP.rileagreeFE1 = feglm(
        appointed ~
          member_pct +
          oppoutin*delta_lr*billdisagree_wo
        |
          1
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.rileagreeFE1)

      mRAP.rileagreeFEcabparty = feglm(
        appointed ~
          member_pct +
          oppoutin +
          oppoutin*delta_lr*billdisagree_wo  |
          cabparty
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.rileagreeFEcabparty)

      mRAP.rileagreeFEcabpartycom = feglm(
        appointed ~
          member_pct +
          oppoutin*delta_lr*billdisagree_wo
        |
          cabparty +
          title_short
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.rileagreeFEcabpartycom)

      mRAP.rileagreeFEcabpartycom_term = feglm(
        appointed ~
          member_pct +
          oppoutin*delta_lr*billdisagree_wo
        |
          cabparty^EPterm +
          title_short^EPterm
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.rileagreeFEcabpartycom_term)

      # Generate table A12
      models_rileagree_RAP = list(mRAP.rileagreeFE1,
                                  mRAP.rileagreeFEcabparty,
                                  mRAP.rileagreeFEcabpartycom,
                                  mRAP.rileagreeFEcabpartycom_term)
      FE_add = data.frame(matrix(c("Party $\\times$ Cabinet FE","", "\\checkmark", "\\checkmark", "\\checkmark",
                                   "Committee FE", "", "", "\\checkmark", "\\checkmark",
                                   "Committee $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Party $\\times$ cabinet $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Log. Lik.", round(unlist(lapply(models_rileagree_RAP, logLik)), 3)),
                                 ncol = length(models_rileagree_RAP) + 1, byrow = T))
      attr(FE_add, "position") = 26 + 1:nrow(FE_add)

      regtab = modelsummary(models_rileagree_RAP,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            width = c(0.4, rep(0.6/length(models_rileagree_RAP), length(models_rileagree_RAP))),
                            title = "Binary logistic regression estimates. DV: Rapporteur assignment. Standard errors in parentheses.\\label{tab:regtab_rileagree_RAP}")# c('*' = .05, '**' = .01, "***" = .001))
      regtab
      save_tt(regtab, here("tabA12.txt"), overwrite = T)
      rm(regtab, FE_add, models_rileagree_RAP)

      # Generate Figure A9
      pdat = plot_slopes(mRAP.rileagreeFEcabpartycom_term, variables = "delta_lr", condition = list("billdisagree_wo", "oppoutin"), draw = F)
      pdat$oppoutin = case_match(pdat$oppoutin,
                                 "opposition" ~ "Opposition",
                                 "outsider" ~ "Outsider",
                                 "insider" ~ "Insider")

      ggplot(pdat, aes(x = billdisagree_wo, y = estimate, ymin = conf.low, ymax = conf.high)) +
        geom_rug(inherit.aes = F,
                 data = d_RAP %>% dplyr::select(billdisagree_wo,oppoutin) %>%
                   mutate(oppoutin = case_match(oppoutin,
                                                "opposition" ~ "Opposition",
                                                "outsider" ~ "Outsider",
                                                "insider" ~ "Insider")),
                 mapping = aes(x = billdisagree_wo),
                 alpha = 0.3) +
        facet_wrap(~oppoutin)+
        geom_ribbon(linetype = 2, color = "grey20", fill = "grey80", alpha = 0.5) +
        geom_line() +
        theme_bw() +
        theme(strip.background = element_blank(),
              strip.text = element_text(size = 11)) +
        labs(x = "Bill disagreement",
             y = TeX("Effect of misalignment ($\\Delta_{policy}$) on rapporteur assignment"),
        )
      ggsave(here("figA9.pdf"), width = 6.5, height = 4)

  ## Misalignment with Insider ----
    ### Committee appointment ----
      # Fit models
      mCOM.policyoutsiderFE1 = feglm(
        member ~
          oppoutin + oppoutin:delta_policy_outsider
        |
          1,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyoutsiderFE1)

      mCOM.policyoutsiderFEcabparty = feglm(
        member ~
          oppoutin + oppoutin:delta_policy_outsider
        |
          cabparty,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyoutsiderFEcabparty)


      mCOM.policyoutsiderFEcabpartycom = feglm(
        member ~
          oppoutin + oppoutin:delta_policy_outsider
        |
          cabparty + title_short,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyoutsiderFEcabpartycom)

      mCOM.policyoutsiderFEcabpartycom_term = feglm(
        member ~
          oppoutin + oppoutin:delta_policy_outsider
        |
          cabparty^EPterm + title_short^EPterm,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_COM
      )
      summary(mCOM.policyoutsiderFEcabpartycom_term)

      # Generate table A13
      models_policyoutsider_COM = list(mCOM.policyoutsiderFE1,
                                       mCOM.policyoutsiderFEcabparty,
                                       mCOM.policyoutsiderFEcabpartycom,
                                       mCOM.policyoutsiderFEcabpartycom_term)
      FE_add = data.frame(matrix(c("Party $\\times$ Cabinet FE","", "\\checkmark", "\\checkmark", "\\checkmark",
                                   "Committee FE", "", "", "\\checkmark", "\\checkmark",
                                   "Committee $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Party $\\times$ cabinet $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Log. Lik.", round(unlist(lapply(models_policyoutsider_COM, logLik)), 3)),
                                 ncol = length(models_policyoutsider_COM) + 1, byrow = T))
      attr(FE_add, "position") = 10 + 1:nrow(FE_add)

      regtab = modelsummary(models_policyoutsider_COM,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            title = "Binary logistic regression estimates. DV: Committee appointment. Standard errors in parentheses.\\label{tab:regtab_deltapolicyoutsider_COM}")# c('*' = .05, '**' = .01, "***" = .001))
      regtab
      save_tt(regtab, here("tabA13.txt"), overwrite = T)
      rm(regtab, models_policyoutsider_COM, FE_add)

      # Generate Figure A10
      pdat = plot_slopes(mCOM.policyoutsiderFEcabpartycom_term, variables = "delta_policy_outsider", by = "oppoutin", draw = F)
      pdat
      pdat$oppoutin = case_match(pdat$oppoutin,
                                 "outsider" ~ "Outsider",
                                 "opposition" ~ "Opposition",
                                 "insider" ~ "Insider")
      ggplot(pdat, aes(x = oppoutin, ydist = dist_normal(estimate, std.error))) +
        geom_hline(yintercept = 0, linetype = 2) +
        stat_eye(fill = "grey80", .width = c(.90, .95), alpha = 0.5) +
        stat_eye(fill = NA, .width = c(.90, .95)) +
        geom_point(aes(y = estimate), size = 2.5) +
        theme_bw() +
        labs(x = "Party status", y = TeX("Effect of misalignment ($\\widetilde{\\Delta}_{policy}$) on committee membership"))
      ggsave(here("figA10.pdf"), width = 6.5, height = 4)

    ### Rapporteur assignment ----
      # Fit models
      mRAP.policyoutsideragreeFE1 = feglm(
        appointed ~
          member_pct +
          oppoutin +
          oppoutin:delta_policy_outsider*billdisagree_wo +
          oppoutin:billdisagree_wo
        |
          1
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyoutsideragreeFE1)
      # NOTE: this model formula is equivalent to the following model. however, this version doesn't falsely plot effect slopes for insider for which delta_policy_outsider is constant.
      mRAP.sanity = feglm(
        appointed ~
          member_pct +
          oppoutin*delta_policy_outsider*billdisagree_wo
        |
          1
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.sanity)

      mRAP.policyoutsideragreeFEcabparty = feglm(
        appointed ~
          member_pct +
          oppoutin +
          oppoutin:delta_policy_outsider*billdisagree_wo +
          oppoutin:billdisagree_wo
        |
          cabparty
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyoutsideragreeFEcabparty)

      mRAP.policyoutsideragreeFEcabpartycom = feglm(
        appointed ~
          member_pct +
          oppoutin +
          oppoutin:delta_policy_outsider*billdisagree_wo +
          oppoutin:billdisagree_wo
        |
          cabparty +
          title_short
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyoutsideragreeFEcabpartycom)

      mRAP.policyoutsideragreeFEcabpartycom_term = feglm(
        appointed ~
          member_pct +
          oppoutin +
          oppoutin:delta_policy_outsider*billdisagree_wo +
          oppoutin:billdisagree_wo
        |
          cabparty^EPterm +
          title_short^EPterm
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )

      mRAP.policyoutsideragreeFEcabpartycom_term = feglm(
        appointed ~
          member_pct +
          oppoutin +
          oppoutin:delta_policy_outsider*billdisagree_wo +
          oppoutin:billdisagree_wo

        |
          cabparty^EPterm +
          title_short^EPterm
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP,
        fixef.keep_names = T
      )
      summary(mRAP.policyoutsideragreeFEcabpartycom_term)

      # Generate table A14
      models_policyoutsideragree_RAP = list(mRAP.policyoutsideragreeFE1,
                                            mRAP.policyoutsideragreeFEcabparty,
                                            mRAP.policyoutsideragreeFEcabpartycom,
                                            mRAP.policyoutsideragreeFEcabpartycom_term)
      FE_add = data.frame(matrix(c("Party $\\times$ Cabinet FE","", "\\checkmark", "\\checkmark", "\\checkmark",
                                   "Committee FE", "", "", "\\checkmark", "\\checkmark",
                                   "Committee $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Party $\\times$ cabinet $\\times$ term FE", "", "", "", "\\checkmark",
                                   "Log. Lik.", round(unlist(lapply(models_policyoutsideragree_RAP, logLik)), 3)),
                                 ncol = length(models_policyoutsideragree_RAP) + 1, byrow = T))
      attr(FE_add, "position") = 22 + 1:nrow(FE_add)

      regtab = modelsummary(models_policyoutsideragree_RAP,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            width = c(0.4, rep(0.6/length(models_policyoutsideragree_RAP), length(models_policyoutsideragree_RAP))),
                            title = "Binary logistic regression estimates. DV: Rapporteur assignment. Standard errors in parentheses.\\label{tab:regtab_policyoutsideragree_RAP}")# c('*' = .05, '**' = .01, "***" = .001))
      regtab
      save_tt(regtab, here("tabA14.txt"), overwrite = T)
      rm(regtab, models_policyoutsideragree_RAP, FE_add)

      # Generate figure A11
      pdat = plot_slopes(mRAP.policyoutsideragreeFEcabpartycom_term, variables = "delta_policy_outsider", condition = list("billdisagree_wo", "oppoutin"), draw = F)
      pdat$oppoutin = case_match(pdat$oppoutin,
                                 "opposition" ~ "Opposition",
                                 "outsider" ~ "Outsider",
                                 "insider" ~ "Insider")
      pdat %>% filter(billdisagree_wo == 1)

      ggplot(pdat, aes(x = billdisagree_wo, y = estimate, ymin = conf.low, ymax = conf.high)) +
        geom_rug(inherit.aes = F,
                 data = d_RAP %>% dplyr::select(billdisagree_wo,oppoutin) %>%
                   mutate(oppoutin = case_match(oppoutin,
                                                "opposition" ~ "Opposition",
                                                "outsider" ~ "Outsider",
                                                "insider" ~ "Insider")),
                 mapping = aes(x = billdisagree_wo),
                 alpha = 0.3) +
        facet_wrap(~oppoutin)+
        geom_ribbon(linetype = 2, color = "grey20", fill = "grey80", alpha = 0.5) +
        geom_line() +
        theme_bw() +
        theme(strip.background = element_blank(),
              strip.text = element_text(size = 11)) +
        labs(x = "Bill disagreement",
             y = TeX("Effect of misalignment ($\\widetilde{\\Delta}_{policy}$) on rapporteur assignment"),
        )
      ggsave(here("figA11.pdf"), width = 6.5, height = 4)

  ## Sensitivity to influential observations ----
      ### Get cook's D ----
      cd = cooks.distance(mRAP.policyagreeFEcabpartycom_term)

      # Generate Figure A12
      ggplot(data.frame(cd), aes(x = cd)) +
        geom_histogram(bins = 100) +
        geom_vline(xintercept = quantile(cd, p = c(0.95, 0.99), na.rm = T)) +
        theme_bw() +
        labs(x = "Cook's distance (Model A4)")
      ggsave(here("figA12.pdf"), width = 6.5, height = 4)

      # Get data for subset models
      d_RAP_cd = d_RAP
      d_RAP_cd$cook = cd

    ### Fit subset models ----
      mRAP.policyagreeFEcabpartycom_term_lowCD99 = feglm(
        appointed ~
          member_pct +
          oppoutin*delta_policy*billdisagree_wo
        |
          cabparty^EPterm +
          title_short^EPterm
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP_cd %>% filter(cook < quantile(cook, p = 0.99, na.rm = T)),
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFEcabpartycom_term_lowCD99)

      mRAP.policyagreeFEcabpartycom_term_lowCD975 = feglm(
        appointed ~
          member_pct +
          oppoutin*delta_policy*billdisagree_wo
        |
          cabparty^EPterm +
          title_short^EPterm
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP_cd %>% filter(cook < quantile(cook, p = 0.975, na.rm = T)),
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFEcabpartycom_term_lowCD975)

      mRAP.policyagreeFEcabpartycom_term_lowCD95 = feglm(
        appointed ~
          member_pct +
          oppoutin*delta_policy*billdisagree_wo
        |
          cabparty^EPterm +
          title_short^EPterm
        ,
        family = binomial,
        vcov = "iid",
        fixef.rm = "none",
        data = d_RAP_cd %>% filter(cook < quantile(cook, p = 0.95, na.rm = T)),
        fixef.keep_names = T
      )
      summary(mRAP.policyagreeFEcabpartycom_term_lowCD95)

      # Generate table A15
      models_policyagree_RAP_CD = list(mRAP.policyagreeFEcabpartycom_term,
                                       mRAP.policyagreeFEcabpartycom_term_lowCD99,
                                       mRAP.policyagreeFEcabpartycom_term_lowCD975,
                                       mRAP.policyagreeFEcabpartycom_term_lowCD95)
      names(models_policyagree_RAP_CD) = c("A4", "A4 [99% low CD]", "A4 [97.5% low CD]",
                                           "A4 [95% low CD]")

      FE_add = data.frame(matrix(c("Party $\\times$ Cabinet FE", "\\checkmark", "\\checkmark","\\checkmark", "\\checkmark",
                                   "Committee FE", "\\checkmark", "\\checkmark","\\checkmark", "\\checkmark",
                                   "Committee $\\times$ term FE", "\\checkmark","\\checkmark", "\\checkmark", "\\checkmark",
                                   "Party $\\times$ cabinet $\\times$ term FE", "\\checkmark", "\\checkmark","\\checkmark", "\\checkmark",
                                   "Log. Lik.", round(unlist(lapply(models_policyagree_RAP_CD, logLik)), 3)),
                                 ncol = length(models_policyagree_RAP_CD) + 1, byrow = T))
      attr(FE_add, "position") = 24 + 1:nrow(FE_add)

      regtab = modelsummary(models_policyagree_RAP_CD,
                            #output = "latex",
                            coef_map = cmap,
                            gof_omit = "AIC|BIC|RMSE|F|Std.Errors|R",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            width = c(0.4, rep(0.6/length(models_policyagree_RAP_CD), length(models_policyagree_RAP_CD))),
                            title = "Binary logistic regression estimates. DV: Rapporteur assignment. Standard errors in parentheses.\\label{tab:regtab_policyagree_RAP_CD}")
      regtab
      save_tt(regtab, here("tabA15.txt"), overwrite = T)
      rm(regtab, models_policyagree_RAP_CD, FE_add)
